Pattern formation in the dipolar Ising model on a two-dimensional honeycomb lattice 
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We present Monte Carlo simulation results for a two-dimensional Ising model with ferromagnetic 
nearest-neighbor couplings and a competing long-range dipolar interaction on a honeycomb lattice. 
Both structural and thermodynamic properties are very similar to the case of a square lattice, with 
the exception that structures reflect the sixfold rotational symmetry of the underlying honeycomb 
lattice. To deal with the long-range nature of the dipolar interaction we also present a simple method 
of evaluating effective interaction coefficients, which can be regarded as a more straightforward 
alternative to the prevalent Ewald summation techniques. 
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I. INTRODUCTION 

The two-dimensional Ising model with ferromagnetic 
nearest-neighbor interactions and long-range antiferro- 
magnetic interactions is probably the simplest model sys- 
tem for the formation of magnetic domains, for instance 
see T. Garel and S. DoniachP Physical systems for which 
the application of such a model is justified are, for exam- 
ple, ultrathin metal films on metal substrates, provided 
that there is a strong magnetocrystalline anisotropy which 
favours spin alignment perpendicular to the plane of the 
filmP In addition to the nearest-neighbor coupling via ex- 
change interaction, the magnetic dipole-dipole interaction 
then realizes a long-range antiferromagnetic coupling that 
decreases as r -3 . A thin magnetic film may therefore be 
described by the following Hamiltonian: 
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where Sf is the z component of a spin 1/2 operator 
on site i, J > corresponds to the ferromagnetic ex- 
change interaction constant, g defines the strength of 
the dipolar coupling and B z denotes an external mag- 
netic field oriented along the z direction. Without loss of 
generality, we measure distances r in units of the nearest- 
neighbor distance. Since the dipolar interaction is inher- 
ently antiferromagnetic (g < 0) and a purely dipolar Ising 
model (J = 0) has the usual antiferromagnetic ground 
state, an additional antiferromagnetic exchange interac- 
tion would simply increase the transition temperatureP 
Therefore, only the case of a ferromagnetic exchange in- 
teraction (J > 0) is interesting. We also want to assume 
for the remainder of this article that the exchange interac- 
tion is strong enough so that it is the dominant coupling 
between nearest neighbors, therefore J > \g\. 

As metal-on-metal films have many technological appli- 
cations, including for example electronics, data storage, 
and catalysisP extensive studies have been performed to 
investigate the properties of Hamiltonian ([!]). Most of 
the metal- on- metal films happen to be square or triangu- 
lar lattice systems, so these studies naturally have had 



their focus on square^H^ and triangular^ lattices. Two- 
dimensional magnets with other lattices are less obvi ous, 
but recently the honeycomb lattice has been discussecP^ 
in the context of the (111) bilayer of LaNiC>3, which shows 
a very rich magnetic phase diagram. To the best of our 
knowledge there are no publications dealing with Hamilto- 
nian ([!]) on an underlying honeycomb lattice. This article 
presents Monte Carlo simulation results for the thermo- 
dynamic and structural properties of such a system. 

Let us briefly recall the existing results for the square 
lattice (for more extensive reviews, see Refs. l2land[TTj): It 
has been analytically established^ that the ground state 
shows a striped pattern of width h, where h increases 
exponentially with the relative strength J/\g\ of the ex- 
change interaction in the limit J> \g\. This implies that 
even an infinitesimal antiferromagnetic dipolar interaction 
will destroy the spontaneous magnetization of the purely 
ferromagnetic ground state. Intermediate between the 
low-temperature striped phase and the high-temperature 
paramagnetic phase with maximum entropy, a third phase 
is found which shows well defined magnetic domains that 
form mazelike patterns P This phase has been called the 
tetragonal phase due to the predominantly rectangular 
corners of those domains. Numerical calculations of the 
structure factor 
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where A: is a wavevector in reciprocal space, show the four- 
fold rotational symmetry of the tetragonal phaseP The 
transition from the striped phase to the tetragonal phase 
is accompanied by a sharp peak in the specific heat at 
temperatures below the Shottky anomaly shoulder, which 
itself can be associated with the tetragonal-paramagnetic 
transition.^ Larger values of J/\g\ generally increase all 
temperatures and make the striped-tetragonal peak less 
pronounced relative to the Shottky shoulderP More recent 
studies^ have revealed the existence of an intermediate 
nematic phase between the striped phase and the tetrag- 
onal phase that is only stable for very specific values 
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of J/\g\ close to some of the ground state stripe width 
transitions. 

In the present work we examine the domain pattern 
formation on a honeycomb lattice in the framework of 
the dipolar Ising model and compare the obtained re- 
sults to the square lattice case. In particular we want 
to determine whether a change of the underlying lattice 
only modifies transition temperatures or also changes 
the system's behavior qualitatively. For this purpose, we 
present an efficient method to numerically evaluate effec- 
tive interaction coefficients. This method can be regarded 
as a simple alternative to the known Ewald summation 
techniques. 




II. METHODS 

Monte Carlo simulations of spin systems with long- 
range interactions are inherently difficult. The reason 
for that lies in the nontrivial implementation of periodic 
boundary conditions. Also, the 0(N 2 ) number of cou- 
plings, where N denotes the number of spins, make the 
evaluation of energy differences computationally expen- 
sive. In the following we will describe the methods we 
used to obtain our results. 

As the tetragonal and striped phases that occur on a 
square lattice clearly reflect the rotational symmetries 
of the underlying lattice, special care should be taken in 
case of the honeycomb lattice to ensure that the shape 
and boundary conditions of the finite system allow the 
formation of structures with the expected rotational sym- 
metries. Fig. [T] illustrates the three rotational symmetries 
of the honeycomb lattice. 

Finite size effects can be reduced through the use of 
periodic boundary conditions. For systems with long- 
range interactions, periodic boundary conditions can be 
implemented by tiling the entire space with replicas of 
the original finite system. 5 Note that this is equivalent 
to the treatment of an infinite system, where only states 
with certain translational invariances are considered. 

Both the rotational symmetries and the requirement 
of the system to be properly tileable are satisfied if one 
chooses the simulated system's shape to be a regular 
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FIG. 1. (Color online) Rotational symmetries of the honey- 
comb lattice: threefold symmetry for rotations around one of 
the basis atoms (point 1), sixfold rotational symmetry around 
the center of the unit cell (point 2) and twofold symmetry 
around the center of nearest-neighbor connections (point 3). 




FIG. 2. Periodic boundary conditions on a honeycomb lattice. 
The original finite system's shape is approximated by the 
central black hexagon. Its replicas form a triangular Bravais 
lattice, whose basis vectors are marked in black. According to 
our nomenclature this would be a (3, 3) aggregation. 



hexagon. This hexagon unit as well as the hexagons' 
tiling is illustrated in Fig. [2] Let us introduce the name 
"aggregation" for the combination of the original system 
and all of its replicas. Note the equivalence between the 
tiling of the system's replicas to form the aggregation and 
the tiling of the unit cells within the original system. It 
is convenient to use the side length of the original system 
in number of unit cells as a measure for the size n of the 
simulated system. Since the aggregation itself is also a 
regular hexagon, it is then straightforward to measure 
its size m as the length of its sides in units of replicas. 
In this way, the entire size and shape of the aggregation 
is specified by a pair (n, m) of integers. Note that the 
original system as well as the entire aggregation are sixfold 
rotationally symmetric around their center, which implies 
threefold and twofold symmetry. 

Now that we have reduced the infinitely large system to 
a finite system with replicas, we can make use of this new 
periodicity. Ignoring for the moment the coupling to an 
external magnetic field, we rewrite the Hamiltonian |T]) by 
using Sf = S z {ri) = S z (R+fi) where R is the translation 
from the original system to any of the replicas: 
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Here the function A is a nonzero constant for nearest 
neighbors (NN) only, whereas T takes care of the dipolc- 
dipole interaction without the unphysical interaction of a 
spin with itself. Note, however, that a spin does interact 
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with its own copies in the replicated systems at R: 



A(fi,R + fj) 



J if f and R + fj are NN, 
otherwise, 







if fi = R + fj 
otherwise. 
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As the sum over all replicas depends only on the indices i 
and j, but not on the orientation of those spins, it can 
be calculated in advance, which reduces the Hamiltonian 
to the general Ising model Hamiltonian with effective 
interaction coefficients: 

H = - E J if s i s l - BZ E s i > ( 4a ) 
J tf = E [ A (^ n + ? i) + r ( f - n + ■ ( 4b ) 



In addition to the exchange interaction, if Sf and 5| 
(or any of its copies) are nearest neighbors, this effective 
interaction coefficient J?? also includes the dipole-dipole 
interaction of Sf with S? and all copies of £| . 

It is obvious that, once the effective interaction coeffi- 
cients have been calculated, the time required to perform 
a Monte Carlo step for a (n, m) aggregation will only de- 
pend on the number of spins and therefore on the size n of 
the original system. The accuracy of the effective interac- 
tion coefficients depends on the size m of the aggregation 
though, but since the Jf^ only have to be calculated once, 
the aggregation size m can be quite large. Note that an 
increase in the system size n for constant values of m also 
leads to an increase in the accuracy of the effective in- 
teraction coefficients, as it corresponds to larger absolute 
values of the translation vectors R and therefore increases 
the distance at which the sum is truncated. 

Due to the finite size of the aggregation, the dipolar 
part of the effective interaction strength will systemati- 
cally be underestimated. Assuming a (n, m) aggregation, 
one can try to compensate for this by calculating the 
effective interaction JSp of a spin with its own copies for 
a (n,m') aggregation, where m' S> m. As Sf will cer- 
tainly not be its own nearest neighbor, the coefficient Jff 
contains only dipolar interactions. In the limit m! — > oo, 
the difference between Jff for (n, m) and (n, m!) aggre- 
gations will just be the remaining dipole interaction that 
has been neglected with the (n, m) aggregation. One can 
now add this difference to all effective interaction coef- 
ficients to compensate approximately for the systematic 
underestimation. Note that this correction, which has 
been calculated for the interaction Jff with the spin's 
own copies, is added to all other coefficients, including Jfj* 
with i ^ j. It therefore also accounts approximately for 
the interaction with copies of all other spins. This approx- 
imation is justified by the fact that the dipolar potential 
is almost flat for very large r and therefore not sensitive 
to the exact relative positions. Because the actual cal- 
culation with the larger (n, mf) aggregation only has to 



be performed for a single spin and its own copies, the 
computational effort is negligible. 

A significant speedup in the calculation of the effective 
interaction coefficients can be achieved if one manages to 
exploit symmetries to reduce the number of coefficients 
that have to be calculated. As only relative positions mat- 
ter, it is quite obvious that not all of the N 2 coefficients 
in a system with N spins will actually be different. In 
principle, one could also reduce the memory footprint of 
the coefficient table in this way, but one would have to 
evaluate very carefully if the overhead of selecting the 
right coefficient does not slow down the Monte Carlo step 
significantly. 

This direct calculation of the effective interaction coef- 
ficients is, compared to the alternative Ewald summation 
techniques,^ a very simple and straightforward method. 
For this particular system it works very well, since it is a 
two-dimensional system where the long-range interaction 
decreases with r~ 3 . This dependence makes the series 
unproblematically convergent. Using the proposed meth- 
ods, we were able to replicate the known results for the 
square lattice within the statistical uncertainty. All of 
our honeycomb lattice simulations were performed on a 
(24, 24) aggregation, which corresponds to a system of 
N = 3314 spins. A (24, 1000) aggregation was used to 
compensate for the underestimation. This calculation 
resulted in a relative increase in the effective coefficients 
of about 10~ 2 for couplings of spins that are far away 
from each other, and about 10 -6 for nearest neighbors. 

Metropolis dynamics^ with uniformly distributed sin- 
gle spin-flip attempts were used to evolve the system. 
Markov chain correlation was dealt with using binning 
and bootstrapping^ techniques. The energy and temper- 
ature scale is defined by measuring J and fc^T in units 
of -g. 



III. RESULTS 

In order to perform a Monte Carlo simulation over 
the whole temperature range, one first has to determine 
which states have the lowest energies. Initializing the 
equilibration period of a low temperature simulation with 
a state that is typical for high temperatures will not 
correctly equilibrate the system: An almost steepest de- 
scent in energy will most likely trap the system in a local 
minimum from which no physical information can be ex- 
tracted. Therefore we have first performed a simulated 
annealing to relax the system to a low energy state, which 
was then used to initialize the equilibration period of the 
simulations. 

We have found the simulated annealing to result in 
a striped state if the system's temperature is decreased 
linearly from T = 5 to over the course of 10 7 Monte 
Carlo steps. This indicates that the ground state of the 
dipolar Ising model on the honeycomb lattice is indeed 
striped, similar to the state that is plotted in Fig. [3]ja). 
Note that one has to be very careful to choose the system 
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FIG. 3. (Color online) (a) Striped phase for J = 6 at ksT = 
0.4 and (b) hexagonal phase for J = 6 at ksT = 1.25. 



35 

30 

25 

20 

15 

10 

5 






u 





kx 



FIG. 4. (Color online) Structure factor cr(fe) in the hexagonal 
phase at ksT — 1.25 and J = 6. 



size to be compatible with the width of the stripes that 
the system would like to form, in order not to introduce 
an artificial frustration. As the width of the stripes is 
determined by the relative strength of exchange and dipo- 
lar interaction, one can also adjust J to make the stripe 
width compatible with the given system size. We have 
found J = 6 to result in a stripe width that is compatible 
with the system size of n = 24 which was used in all 
our simulations. Naturally, the stripes become wider if 
the strength of the ferromagnetic exchange interaction is 
increased. 

Having determined the low-temperature states to be 
striped by simulated annealing, we have equilibrated 
the system at constant temperatures for 2 x 10 5 Monte 
Carlo steps and recorded thermodynamic and structural 
properties for 8 x 10 6 time steps. We find that the low- 
temperature striped phase is followed by a phase with a 
complex domain structure which is plotted in Fig. |3]jb). 
We call this phase hexagonal, due to the sixfold rotational 
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FIG. 5. (Color online) Specific heat capacity per spin c(ksT) 
as a function of temperature for a relative interaction strength 
of J = 6. The line is just a guide to the eye. 



symmetry that is visible in the structure factor cr(fc) plot- 
ted in Fig. [4] The six-peaked shape of the structure factor 
remains unchanged if the sum in Eq. ([2| is restricted to 
one of the two basis atoms in the honeycomb unit cell. We 
can conclude from this that the hexagonal phase shows 
the rotational symmetry of the triangular Bravais lattice 
that underlies the honeycomb lattice. It is interesting to 
note that a change of the lattice from square to honey- 
comb changes only the intermediate-temperature phase 
from tetragonal to hexagonal, while the low-temperature 
phase is striped for both lattices. 

We do not observe the three-peaked structure of the 
energy histrogram that was used by Cannas et alP to 
identify the square lattice's nematic phase. We therefore 
have to conclude that there is no nematic phase for the 
particular value of J/\g\ used in our calculations. 

Similar to the striped-tetragonal transition on the 
square lattice, the striped-hexagonal transition manifests 
itself in a sharp peak in the specific heat capacity at tem- 
peratures below the expected Shottky anomaly shoulder 
as shown in Fig. [5j We expect the striped-hexagonal peak 
to become less pronounced for larger J, as is the case on 
a square lattice. 

In summary, by means of Monte Carlo simulations we 
have investigated the temperature dependence of domain 
pattern formation in a dipolar Ising model on a two- 
dimensional honeycomb lattice. Taking into account the 
symmetries of the honeycomb lattice, we have reduced 
the infinite system to a finite system with replicas, which 
can be described by the general Ising model Hamiltonian 
with effective interaction coefficients. We have presented 
a straightforward method of calculating the effective in- 
teraction coefficients and a procedure to compensate ap- 
proximately for their systematic underestimation. We 
find that the honeycomb lattice shows two distinct phase 
transitions: from a striped phase at low temperatures 
via a hexagonal phase at intermediate temperatures to a 
disordered phase at high temperatures. Both transitions 
are associated with maxima in the specific heat. While 
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the thermodynamic properties of the honeycomb lattice 
system are found to be virtually identical to its square lat- 
tice counterpart, the emerging patterns clearly reflect the 
different rotational symmetry of the underlying lattice. 

Future work should focus on obtaining the entire phase 
diagram of the model in order to determine the values 
of J/\g\ at which the ground state shows transitions in 



stripe width. For those values one can then systematically 
search for the honeycomb lattice's analogon to the nematic 
phase^l in between the striped and the hexagonal phases. 
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